Screw dislocation in zirconium: an ab initio study 
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Plasticity in zirconium is controlled by 1/3(1210) screw dislocations gliding in the prism planes 
of the hexagonal close-packed structure. This prismatic and not basal glide is observed for a given 
set of transition metals like zirconium and is known to be related to the number of valence electrons 
in the d band. We use ab initio calculations based on the density functional theory to study the 
core structure of screw dislocations in zirconium. Dislocations are found to dissociate in the prism 
plane in two partial dislocations, each with a pure screw character. Ab initio calculations also show 
that the dissociation in the basal plane is unstable. We calculate then the Peierls barrier for a 
screw dislocation gliding in the prism plane and obtain a small barrier. The Peierls stress deduced 
from this barrier is lower than 21 MPa, which is in agreement with experimental data. The ability 
of an empirical potential relying on the embedded atom method (EAM) to model dislocations in 
zirconium is also tested against these ab initio calculations. 

PACS numbers: 61.72.Lk, 61.72.Bb 
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I. INTRODUCTION 

Plasticity in a-zirconium is controlled by dislocations 
with a 1/3(1210) Burgers vector gliding in the prism 
planes of the hexagonal compact (hep) latticed— The rel- 
ative ease of prismatic glide compared to basal glide has 
been shown to be linked to the ratio of the correspond- 
ing stacking fault energies, which in turn is controlled by 
the electronic structured In particular, Legrand 5 used a 
tight binding model to show that prismatic slip in tran- 
sition metals of the IV B column (Zr, Ti, Hf ) originates 
from the electronic filling of the valence d band. As a 
consequence, it appears necessary to take into account 
the anisotropy of the d orbital, and hence the angular 
dependence of the atomic bonding, to model dislocations 
in these transition metals^. One cannot therefore rely on 
central forces empirical potential and needs to take ac- 
count of the electronic structure. Tight binding models^ 
or ab initio calculations 9-13 show indeed that a 1/3(1210) 
screw dislocation, either in Zr or in Ti, spreads in prism 
planes, in agreement with the prismatic glide observed 
experimentally. But none of these atomistic simulations 
calculate the Peierls stress of a screw dislocation. Ac- 
cording to some authors j 7 -^ its core structure is not com- 
pletely planar, which may be the cause of a high Peierls 
stress. 

It is true, experimentally, that screw dislocations 
glide with difficulty compared to other dislocation char- 
acters in zirconium or titanium alloys: characteris- 
tic microstructures, with long and straight dislocations 
aligned along their screw orientation, are observed at 
low temperature^i^iii^— and the flow stress is strongly 
temperature dependent 1 6 1 1 9 — in agreement with the 
assumption of a high Peierls barrier which must be over- 
come by the nucleation of double kinks. But experiments 
also show that the yield stress in zirconium or in tita- 
nium strongly decreases with a decreasing amount of in- 
terstitial impurities like oxygeni 2 i 3 i 16 ' 19 i 22 : 23 It is there- 
fore probable that the high Peierls stress of screw dis- 



locations has an extrinsic cause. In pure zirconium or 
pure titanium, this Peierls stress may not be as high and 
screw dislocations are probably gliding as easily as other 
orientations. 

Recently, Mendelev and Ackland2i developed an em- 
pirical potential for Zr in the Embedded Atom Method 
(EAM) formalism. Using this potential, Khater and 
Bacon 2 - showed that it leads to a screw dislocation that 
spontaneously spreads in the prism plane, the configu- 
ration dissociated in the basal plane being metastable. 
They also showed that the Peierls stress of a screw dislo- 
cation gliding in the prism plane is not so different from 
the Peierls stress of an edge dislocation and that this 
stress is small (22 MPa for the screw and 16 MPa for the 
edge). These results therefore support experimental find- 
ings stating that screw dislocations are gliding in prism 
planes with a low Peierls stress in pure zirconium. But 
these simulations rely on a central forces potential, which 
is not well suited to describe dislocations in a hep tran- 
sition metal like Zr, as described above. More reliable 
atomistic simulations, incorporating a description of the 
electronic structure, are therefore needed to confirm this 
easy glide of screw dislocations in pure zirconium. 

This article aims to use ab initio calculations so as 
to fully characterize the core structure of a 1/3(1210) 
screw dislocation in zirconium and estimate its Peierls 
stress. We also examine generalized stacking faults as 
dislocation core structures are closely related to them. 
Two different ab initio methods, Siesta 2 ^ and Pwscp 2 ^, 
have been used. Siesta offers the advantage of effi- 
ciency, allowing simulating more atoms than with a stan- 
dard ab initio code, whereas Pwscf offers the advantage 
of robustness. All calculations are also performed with 
Mendelev and Ackland EAM potential^ so as to iden- 
tify its ability to model dislocations in zirconium. In 
addition, this empirical potential is used to study the 
convergence of our results with the size of the simulation 
cell. 
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II. ATOMIC INTERACTION MODELING 

Atomistic simulations have been performed both with 
an empirical interatomic potential and with ab initio cal- 
culations. The empirical potential that we used is the 
EAM potential developed by Mendelev and Acklandi^ 
This potential is labeled #3 in Ref. HJ. It is supposed to 
be well suited to model dislocations, as ab initio values^ 
of the stacking fault energies in the basal and prism 
planes have been included in the fitting procedure. Us- 
ing this potential, Khater and Bacor*^ showed that a 
1/3 (1210) screw dislocation spontaneously dissociates in 
the prism plane and that a metastable configuration dis- 
sociated in the basal plane also exists. 

The ab initio calculations are relying on the Density 
Functional Theory (DFT) in the Generalized Gradient 
Approximation (GGA) with the functional proposed by 
Perdew, Burke, and Ernzerhof (PBE) and the pseudopo- 
tential approximation. Two different ab initio codes are 
used, Siesta^ and PwsCFi^l 

In the Siesta code^ valence electrons are described 
by a localized basis set corresponding to a linear com- 
bination of pseudoatomic orbitals with 13 functions per 
atom. We used a norm conserving pseudopotential of 
Troulliers-Martins type with 4p electrons included as 
semicore. Electronic density of state is broaden with the 
Mcthfcssel-Paxton function with a broadening of 0.3 eV 
and the integration is performed on a regular grid of 
14 x 14 x 8 k-points for the primitive hep unit cell and an 
equivalent density of k-points for the supercells used for 
the defect calculations. The charge density is represented 
on a real space grid with a spacing of 0.08 A (Mesh cutoff: 
450 Ry) that is reduced to 0.06 A (800 Ry) for dislocation 
calculations. This approach, i.e. the basis and the pseu- 
dopotential, has already been used to study vacancy dif- 
fusion in zirconium^ and comparison with plane waves 
DFT calculations has led to a reasonable agreement. 

In the Pwscf code^i valence electrons are described 
with plane waves using a cutoff energy of 28 Ry. The 
pseudopotential is ultrasoft of Vanderbilt type with 4s 
and 4p electrons included as semicore^ The same k- 
point grid and the same electronic broadening are used 
as with Siesta code. 

To validate these different atomic interaction models, 
it is worth comparing their results to available experi- 
mental data for some bulk properties of Zr. All models 
lead to an equilibrium lattice parameter and a c/a ratio 
in good agreement with experimental data (Tab. In 
particular, a ratio lower than the ideal y/8/3 ~ 1.633 
value is obtained in all cases. 

We also compared the theoretical elastic constants 
with experimental data (Tab.HJ: a good agreement is also 
obtained. The computed elastic constants are the relaxed 
ones^ as the hep lattice contains two atoms in its primi- 
tive unit cell, some internal degrees of freedom may exist 
when applying a homogeneous strain. One needs to allow 
for atomic relaxations when computing Cn, C\2 or Cqq 
constants^ It is also possible to calculate inner elastic 



TABLE I. Bulk properties of hep Zr calculated with differ- 
ent atomic interaction models and compared to experimental 
data: lattice parameter a, c/a ratio of the hexagonal lattice, 
relaxed elastic constants dj , inner elastic constants^ etj and 
dij, phonon frequencies oji and 013 of the optical branches at 
the r point (Eq. [T}, inner elasticity contribution to elastic 
constant 8C12 (Eq. [2J). 





Expt. 


EAM 


Siesta 


PWSCF 


a (A) 


3.232' 12 


3.234 


3.237 


3.230 


c/a 


1.603 32 


1.598 


1.613 


1.601 


Cn (GPa) 


155.4^ 


142. 


140. 


140. 


C 33 (GPa) 

00 \ J 


172.5^ 


168. 


168. 


168. 


C12 (GPa) 


67.2^ 


75. 


86. 


70. 


C 13 (GPa) 


64.6^ 


76. 


68. 


65. 


C44 (GPa) 


36. 3^ 


44. 


24. 


26. 


C 66 (GPa) 


44. 1* 


33.5 


27. 


35. 


en (meV A" 5 ) 




27.0 


17.0 


18.6 


e 33 (meV A" 5 ) 




118. 


122. 


101. 


tfei (meV A" 4 ) 




30.0 


38.9 


36.6 


cji (THz) 


2.66 ±0.0222 


2.60 


2.08 


2.16 


0J3 (THz) 


4.23±0.15 3s 


5.43 


5.57 


5.03 


SC 12 (GPa) 




5.33 


14.3 


11.5 



a Experimental elastic constants 3 ^ have been measured at 4K. 



constants to characterize these internal degrees of free- 
dom. These are also given in Tab. U using the notations 
introduced by Cousins.— Two of these inner elastic con- 
stants, en and 633, are related to the phonon frequencies 
of the optical branches at the T poinli^i 

Wi = 2,/^i, (l) 

V m 

where VL = a 2 c\fS/A is the atomic volume and m the 
atomic mass. The last inner elastic constants (L21 cou- 
ples the internal degrees of freedom to the homogeneous 
strain. It leads to a contribution to the elastic constants 31 

SC 12 = (2) 
en 

If C% are the unrelaxed elastic constants, i.e. the elastic 
constants calculated by imposing a homogeneous strain 
to the hep lattice without letting atoms relax from their 
initial positions, the true elastic constants are given by 31 - 
C\x = C u — dC'12, C12 = Ci 2 + SC12, and Cqq = Cg 6 — 
<5Ci2, all other elastic constants being unchanged. 

III. STACKING FAULT ENERGIES 

Dislocation dissociation is controlled by the existence 
of a metastable stacking fault for the corresponding 
plane. According to the results obtained by Khater and 
Bacon2£ with Mendelev and Ackland EAM potential^ 
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TABLE II. Stacking fault energies in the basal plane, 7b, and 
in the prism plane, y p , calculated with different atomic inter- 
action models, including Vasp calculations of Domain et aL— 
and of Udagawa et aL— 7? is the ratio defined by Legrand, 5 
and d b q and dp q are the dissociation lengths of a screw dislo- 
cation respectively in the basal (Eq. [5} and in the prism plane 
(Eq.©. 



EAM 



Siesta Pwscf Vasf 



Vasi 



7b (mJm~ 2 ) 


198. 


199. 


213. 


200. 


227 


7 P (mJm -2 ) 


135. (274.£ 


233. 


211. 


145. 


197 


R = C 6 67b/C447p 


1.12 


0.96 


1.36 


1.85 


2.1 


< q (A) 


4.0 


2.0 


2.7 


3.4 


3.2 


rfp q (A) 


7.8 


4.6 


5.9 


9.6 


7.4 


a For the EAM calculation of the prismatic 


stacking 


fault 


energy, 



the value in parenthesis corresponds to the maximum in 
a/6[1210], whereas the lower value corresponds to the true 
minimum in a/6 [1210] + 0.14c [0001]. 



a 1/3 (1210) screw dislocation can dissociate either in a 
basal or in a prism plane. To characterize these even- 
tual dissociations, we compute generalized stacking fault 
energie o 14 i 36 - or 7-surfaces - for both the basal and the 
prism planes. 



(at Pwscf - basal plane 
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(c) EAM - basal plane 
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A. Methodology 



7-surfaces describe the energy variation when two 
parts of a crystal are rigidly shifted for different fault 
vectors lying in a given crystallographic plane. Atoms 
are allowed to relax in the direction perpendicular to the 
fault plane. We calculate these 7-surfaces for both the 
basal and prism planes using full periodic boundary con- 
ditions. To introduce only one fault in the simulation 
cell, the same shift as the one corresponding to the fault 
vector is applied to the periodicity vector perpendicular 
to the fault plane. No free surfaces is therefore intro- 
duced in the simulation cell, which allows a fast conver- 
gence of the fault energies with the number of stacked 
planes. A periodic stacking of at least 16 {0001} planes 
is used for the basal fault and 12 {1010} planes for the 
prismatic fault. This corresponds to a distance between 
fault planes /loooi = 8c ~ 41 A and h 1Q i Q = 6\/3a ~ 34 A. 
Increasing the number of planes in the stacking modifies 
the energies by less than 1 mJ m -2 . Generalized stacking 
fault energies are calculated on a regular grid of 10 x 10 
fault vectors and are then interpolated with Fourier se- 
ries. 




1010] 



a/3 [1210] 



FIG. 1. Generalized stacking fault energy in the basal plane 
calculated with (a) Pwscf, (b) Siesta and (c) EAM poten- 
tial. The arrows indicate Burgers vectors of the partial dis- 
locations corresponding to a dissociation in the basal plane. 
The dashed line is the [1100] direction used in Fig. [2] Contour 
lines are drawn at the base every 50mJm -2 . 



B. 7-surfaces 

1. Basal plane 

7-surfaces corresponding to the basal plane are shown 
in Fig. [1] for the three interaction models we used. In all 
cases, a local minimum can be found at 1/3[1100] which 
corresponds to the intrinsic I2 fault. 38 This minimum 
does not vary when full atomic relaxations are allowed 
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a [1100] 

FIG. 2. Generalized staking fault energy in the basal plane 
along the [1100] direction (cf. Fig.[TJ calculated with Pwscf, 
Siesta and EAM potential. 



instead of being constrained to the direction perpendicu- 
lar to the fault plane. All methods give a close value for 
the fault energy 7b in this minimum (Tab. |TTJ) . A good 
agreement is also obtained with the values calculated by 
Domain et al£2- and by Udagawa et alM- using Vasp ab 
initio code. The depth of this local minimum is more 
pronounced with the empirical EAM potential [Fig.QJc)] 
than in ab initio calculations [Fig.flja) and (b)]. This ap- 
pears clearly in Fig. [5] where the fault energy predicted by 
the different interaction models are compared along the 
[1100] direction. We will see latter that this has conse- 
quences on the stability of a screw dislocation dissociated 
in the basal plane. 

The 7-surface calculated with the EAM empirical po- 
tential has another minimum located in 2/3 [1100]. This 
is an artifact of the potential: one expects instead a 
maximum as this fault vector transforms the original 
BABABA stacking of the basal planes in a BABBCB 
stacking. Ab initio calculations confirm that this fault 
vector gives a maximum. Finally, it is worth pointing 
that both ab initio techniques give a very similar 7- 
surfacc: the shapes are identical and the amplitudes do 
not differ by more than 10%. 



2. Prism plane 

7-surfaces for the prism plane are shown in Fig. [3] 
Both ab initio methods lead to a similar 7-surface 
[Figs. G^a) and (b)]. In particular, both Pwscf and 
Siesta predicts the existence of a minimum at halfway of 
the Burgers vector, i.e. in 1/6 [1210]. Like for the basal 
fault, this minimum does not vary when full atomic re- 
laxations are allowed. As can be seen on the projection 
of these 7-surface along the [1210] direction (Fig. 2]), this 
minimum is a little more pronounced with Pwscf than 
with Siesta. The same minimum was also present in the 
Vasp calculations of Domain et a/.,— but they obtained a 
lower value 7 P of the fault energy in this point (Tab. [[[]). 



Ia| Pwscf - prism plane 




c [0001] 



(h) Siesta - prism plane 




c [0001] 



ai$ [1210] 



(e) EAM - prism plane 




c [0001] 



a/3 [1210] 



FIG. 3. Generalized stacking fault energy in the prism plane 
calculated with (a) Pwscf, (b) Siesta and (c) EAM poten- 
tial. The arrows indicate Burgers vectors of the partial dis- 
locations corresponding to a dissociation in the prism plane. 
Contour lines are drawn at the base every 75mJm~ 2 . 



On the other hand, Udagawa et alM- obtained a value 
close to our result using also Vasp ab initio code with 
the PAW method and the PBE-GGA functional. They 
pointed out that the discrepancy arises from an insuf- 
ficient number of stacked planes in the 7-surface calcu- 
lation of Domain et al. The energy of the metastablc 
stacking fault energy in the prism plane appears there- 
fore higher than the value 145 mJm" 2 initially suggested 
by Domain et al.£& both our Pwscf 3 ^ and Siesta cal- 
culations, as well as the Udagawa et al. result, 37 leads to 
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FIG. 4. Generalized staking fault energy in the prism plane 
along the [1210] direction calculated with Pwscf, Siesta and 
EAM potential. 



a value of about 200 mJm -2 . 

The 7-surface calculated with the EAM potential 
is quite different [Fig. [3[c)] as the point located in 
1/6 [1210] is indeed a maximum and not a minimum like 
with ab initio calculations. A minimum is found for a 
point located in a/6 [1210] + 0.14c [0001]. One therefore 
expects that a 1/3(1210) {1010} dislocation dissociates 
into two partial dislocations with a Burgers vector com- 
ponent orthogonal to the one of the perfect dislocation. 
In particular, a screw dislocation should dissociate in two 
partial dislocations with a small edge character. Khater 
and Bacon 25 showed that the empirical potential of Ack- 
land et al— suffers from the same artefact. As noted by 
Bacon and Vitekj^I all central forces potentials stabilize 
indeed such a stacking fault a/6 [1210] + ac [0001] with 
a ^ 0, a minimum also predicted by a simple hard sphere 
model.— This minimum either disappears or is located 
exactly in a/6 [1210] (a = 0) only when the angular de- 
pendence of the atomic interactions is taken into account. 
The value 7 P of the stacking fault energy obtained with 
Mendelev and Ackland EAM potential 24 for this mini- 
mum is much lower than our ab initio value (Tab. HI]) , 
This is quite normal as Mendelev and Ackland used the 
ab initio value of Domain et al. to fit their potential. 



C. Dislocation dissociation 

Before using atomistic simulations to obtain disloca- 
tion core structures and their associated Peierls stress, it 
is worth looking what can be learned from these 7-surface 
calculations. Legrand 5 proposed a criterion based on the 
elastic constants and the metastable stacking fault en- 
ergies to determine if glide occurs in the base or prism 
plane in a hep crystal. According to this criterion, pris- 
matic glide is favored if the ratio R — C667b/C447p is 
larger than 1. Tab \H\ shows that this is the case for the 
EAM potential and the Pwscf calculation, as well as for 
the Vasp calculation of Domain et al£& and of Udagawa 



et al.M- On the other hand, Siesta leads to a value too 
close to 1 to be able to decide between basal and pris- 
matic glide. 

One can also use dislocation elasticity theory 4 ^ to com- 
pute the dissociation distance of a dislocation both in the 
basal and prism planes. According to elasticity theory, 
the energy variation caused by a dissociation of length d 
is 



(3) 



where bW and b^ 2 ) are the Burgers vectors of each par- 
tial dislocation, 7 the corresponding stacking fault en- 
ergy, K the Stroh matri x 44 ' 45 controlling dislocation elas- 
tic energy, and r c the core radius. Minimization of this 
energy leads to the equilibrium dissociation length 



d eq = 



1 



(4) 



When the hep crystal is oriented with the x, y, and z 
axis respectively along the [1010], [0001], and [1210] di- 
rections, for a dislocation lying along the z direction, the 
K matrix is diagonal with its components given by 4 ^ 4849 



1 / - \ C44 (Cn — C13J 

KU = to (Cn + Cl3) V C 33 (C 11+ C 13 + 2C 4i ) 



^22 = \l -^Ku, 



K33 — 2~~; y 2^ 44 ^ n — ^ 12 )' 



where C n = VCiiC 3 3- 

The 7-surface of the basal plane (Fig. [T]) indicates a 
possible dissociation 1/3[1210] -> 1/3[1100] + 1/3[0110]. 
The dissociation length in the basal plane is then, for a 
1/3 [1210] screw dislocation, 



(3X 33 - K n ) 
12 7 b 



(5) 



According to the minimum of the 7-surface in the 
prism plane (Fig. [3]), a 1/3[1210] dislocation can dissoci- 
ate in this plane in two partial dislocations with Burgers 
vectors 1/6[1210] ± ac/a[0001]. The parameter a con- 
trols the position of the stacking fault minimum along 
the [0001] direction, i.e. a — for Pwscf and Siesta, 
and a = 0.14 for the EAM potential. The dissociation 
length of a screw dislocation in the prism plane is then 



d cq = 

Up — 



(K 33 a 2 - 4a 2 if 22 c 2 
47 P 



(0) 



The dissociation lengths a^ q and dp q calculated from 
the elastic constants and the stacking fault energies are 
given in table [TTJ Elastic constants predicted by the 
atomic interaction models are used in each case. For all 
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energy models, one expects a larger dissociation in the 
prism than in the basal plane. We will compare in the 
following section these dissociation lengths predicted by 
elasticity theory with the ones observed in our atomistic 
simulations of the dislocation core structure. 



IV. SCREW DISLOCATION CORE 
A. Methodology 



a) O arrangement b) S arrangement 




® *~ e x || [10101 

e z || [1210] 



dipole is obtained by a tt/2 rotation of D. 

The dislocation dipole is introduced in the simulation 
cells by applying to all atoms the elastic displacement 
predicted by anisotropic elasticity theor y 44 i 45 taking full 
account of the periodic boundary conditions^ A homo- 
geneous strain is also applied to the simulation cell so as 
to cancel the plastic strain introduced by the dislocation 
dipole and minimize the elastic energy. This strain is 
given by 

o _ bAj + bjAj 
ij 2S 

where S = \(Ui A U%) • e z \ is the surface of the simulation 
cell perpendicular to the dislocation lines. This homoge- 
neous strain adds some tilt components to the periodicity 
vectors. Atoms are then relaxed until all components of 
the atomic forces are smaller than 5 meV A -1 for Siesta, 
2 meV A" 1 for Pwscf, and 0.1 meV A" 1 for the EAM 
potential. 

B. Core structure 



FIG. 5. Screw dislocation periodic arrangements used for 
atomistic simulations. Ui and U2 are the periodicity vectors 
of the arrangement, and A the cut vector of the dislocation 
dipole. The thin vertical line corresponds to the prism glide 
plane. 

Our atomistic simulations of the core structure of 
a screw dislocation are based on full periodic bound- 
ary conditions. 50 This requires introducing a dislocation 
dipole in the simulation cell. Two periodic arrangements 
of the dislocations have been used (Fig. [5|). In the O 
arrangement, dislocations with opposite Burgers vectors 
are located on the same prism and basal planes, i.e. the 
two foreseen glide planes. On the other hand, only dislo- 
cation with the same Burgers vectors can be found on a 
given prism or basal plane in the S arrangement. 

Periodicity vectors of the O arrangement are, before 



(a) Pwscf 



(b) Siesta 



(c) EAM 



introducing the dislocations, U\ 



[1010]-mc[0001], 



U 2 = na\ [1010] + mc [0001], and U 3 = a± [1210]. The 
integers n and m are taken equal to keep an aspect ratio 
close to a square. This simulation cell has been used for 
ab initio calculations with n = 4 (128 atoms), n — 5 (160 
atoms), and n — 6 (288 atoms). 



For the S 



arrangement, 
1 



Ux 



na-- 



[10101, u 2 



"2 L-"-"Ji ^2 - 

mc [0001], and U3 — a-| [1210]. Ab initio calculations 



have been performed with n — m and varying between 
n = 5 (100 atoms) and n — 8 (256 atoms). 

Both dislocation arrays are quadrupolar: the vec- 
tor joining the two dislocations composing the primitive 
dipole is D = l/2(Ux + Ui)- Because of the centrosym- 
metry of this arrangement and the symmetry properties 
of the Volterra elastic field, this ensures that the stress 
created by other dislocations is minimal at each disloca- 
tion position. The cut vector A defining the dislocation 



• o • 




[0001] 



jioio] . \ 




®[1210] • 



FIG. 6. Differential displacement maps around one of the two 
1/3[1210] screw dislocations composing the dipole for the S 
periodic arrangement with n = m = 7 (196 atoms). Atoms 
are sketched by circles with a color depending on the (1210) 
plane to which they belong. The arrow between two atomic 
columns is proportional to the [1210] component of the differ- 
ential displacement between the two atoms. Displacements 
smaller than 0.1& are not shown. Crosses x correspond to 
the positions of the two partial dislocations, and + to their 
middle, i.e. the position of the total dislocation. 

Starting from perfect dislocations, atom relaxation 
leads to dislocations spread in the prism plane, whatever 
the interaction model (EAM, Siesta, or Pwscf) and 
whatever the simulation cell used. This can be clearly 
seen by plotting differential displacement maps as intro- 
duced by Vitek. 52 These maps (Fig. |6]) show that the 
strain created by the screw dislocation spreads out in 
the (1010) prism plane and that displacements outside 
this plane are much smaller. 

To characterize this spreading in the (1010) prism 
plane, we extract from our atomistic simulations the dis- 
registry D(x) created by the dislocation. This is de- 
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-10 10 

Position x in prismatic plane (A) 

FIG. 7. Disregistry D(x) created by the screw dislocation 
in the prism plane and corresponding dislocation density 
p(x) = dD(x)/dx for the S periodic arrangement with 196 
atoms (n = m = 7). Symbols correspond to atomistic sim- 
ulations and lines to the fit of the Peierls-Nabarro model to 
these data. For clarity, disregistries D(x) have been shifted 
by 0.25 between the different interaction models. 



fined as the displacement difference between the atoms 
in the plane just above and those just below the dis- 
location glide plane. The derivative of this function, 
p(x) — dD/dx corresponds to the dislocation density. 
Fig. [7] shows the disregistry obtained for the three inter- 
action models. In all cases, the b discontinuity created 
by the screw dislocation does not show a sharp interface, 
but spreads on a distance ~ 10 A. 

Peierls^ and Nabarro^ built a model that leads to 
a simple expression of the disregistry. The analytical 
expression they obtained^ can be extended to a dissoci- 
ated dislocation. As suggested by the prismatic 7-surfacc 
(Fig. |3]), we assume that the screw dislocation dissociates 
in two equivalent partial dislocation separated by a dis- 
tance d. Based on the Peierls-Nabarro model, we write 
the disregistry created by a single dislocation 



D A i s \ (x) =— < archil) 



x — x — d/2 



arctan 



x — xq + d/2 



c 



tt/2 



(7) 



period L is 

00 

D L {x)= ^ D dis i Q (x-nL) 



n— — 00 



= — < arctan 
2vr 



tan(f [x -x - d/2}) 



tanhf^ 



x — xo — d/2 1 
L + 2 



arctan 



tan(f [x-Xo+d/2]) 



tanh 

x — xq + d/2 1 



(8) 



where [_ - J is the floor function. For the O arrangement 
(Fig. [5£i), the disregistry in the prism plane should be 
given by D{x) = Dl(x) — Dl(x — L/2) with L = 2mc, 
whereas it should be D(x) — Dl(x) with L = mc for the 
S arrangement (Fig. \Sjp). 

We fit this analytical expression of the dislocation dis- 
registry to the data coming from the atomistic simula- 
tions. Fig. [7] shows a good agreement between atomistic 
simulations and the model, using only three fitting pa- 
rameters: the dislocation position xq, the dissociation 
length d and the spreading (. This procedure therefore 
allows us to determine the location of the dislocation cen- 
ter. For all interaction models, we find that this center 
lies in between two (0001) atomic planes. One can see in 
Fig.[6]that this position corresponds to a local symmetry 
axis of the differential displacement map. This is differ- 
ent from the result obtained by Ghazisaeidi and Trinkle 13 
in Ti where the center of the screw dislocation was found 
to lie exactly in one (0001) atomic plane, a position that 
corresponds in Zr to the saddle point between two Peierls 
valleys, as it will be shown below. 

The dissociation length of the screw dislocation ob- 
tained through this fitting procedures are shown in Fig. [8] 
for both periodic arrangements. We observe variations 
with the dislocation periodic arrangement used in the 
simulation, as well as with the size of the simulation 
cell. The results are nevertheless close to the predictions 
of elasticity theory f ^III CI) for all the three interaction 
models. We could even see with the EAM potential that 
the dissociation lengths extracted from atomistic simu- 
lations converge to the value given by elasticity theory 
for large enough unit cells. It is thus relevant to describe 
the screw dislocation as dissociated in two partial disloca- 
tions linked by a stacking fault, although the dissociation 
length remains small. 



where xq is the dislocation position, d its dissociation 
length, and £ the spreading of each partial dislocation. 
We need then to take into account that we do not have 
only one dislocation on a given prism plane but a periodic 
array (Fig. [5]). The disregistry created by an array of 



C. Core energy 

We obtain the dislocation core energy by subtracting 
the elastic energy from the excess energy given by the 
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FIG. 8. Dissociation length of the screw dislocation. Sym- 
bols correspond to data extracted from atomistic simulations 
through the fit of the disregistry (Eq. [8} for both the O and 
S periodic arrangements and for different sizes of the simu- 
lation cell. The solid lines are the predictions of elasticity 
theory based on the stacking faults (Tab. [II]). On the right 
vertical axis, the experimental c lattice parameter has been 
used to normalized the dissociation length by the distance 
Ap = c/2 between Peierls valleys. 
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FIG. 9. Core energy of the screw dislocation obtained for 
different sizes of the atomistic simulation cell and for both 
the O and S periodic arrangements (r c = b). 



atomistic simulations. This excess energy is the energy 
difference between the simulation cell containing the dis- 
location dipole and the same cell without any defect. 
The elastic energy calculation takes into account the 
elastic anisotrop y 44 ' 45 and the effect of periodic bound- 
ary conditions^ The obtained core energy are shown in 
Fig. [SJ Like for the dissociation length, some variations 
with the size of the simulation cell and the periodic ar- 
rangement can be observed. We did not manage to link 
both quantities, i.e. the core energy and the dissociation 
length. They are not simply related by the expression of 
the energy variation with the dissociation length (Eq. [3J , 
and the change of the elastic interaction between disloca- 
tions caused by their dissociation could not fully explain 
the variation of the core energy either. As the spread- 



ing of partial dislocations also depends on the size of the 
simulation cells, the variation of the core energy proba- 
bly have a more complex origin than simply a variation 
of the dissociation length. 

We could use, with the EAM potential, a simulation 
cell large enough to obtain a converged value of the core 
energy, 163meVA _1 for a core radius r c = b. It is worth 
pointing that this value does not depend on the periodic 
arrangement used in the simulation (Fig. [9]). This can be 
achieved thanks to a proper account of the core traction 
contribution to the elastic energy^ A difference of ~ 
lOmeV A -1 would have been observed between the S and 
O periodic arrangement without this contribution. Ab 
initio calculations lead to a dislocation core energy of 
145 ± 5meVA _1 for Siesta and 125 ± 5meVA _1 for 
Pwscf (r c = b in both cases). 



D. Dissociation in the basal plane 

Basal slip is observed experimentally only at high tem- 
perature (above 850 K) and for a higher resolved shear 
stress than the one needed to activate prismatic slip<££ At 
lower temperatures, no basal slip could be observed^^ 
even when the monocrystal was oriented so as to favor 
basal slip. Our atomistic simulations lead to a screw 
dislocation configuration dissociated in the prism plane, 
which clearly cannot glide easily in the basal plane. It is 
worth looking if another configuration, which could glide 
in this basal plane, also exists. 
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FIG. 10. Differential displacement map of the metastable 
configuration of a 1/3[1210] screw dislocation dissociated in 
the basal plane, as obtained with the EAM potential for the 
O periodic arrangement with n = 6 (288 atoms). Crosses x 
correspond to the positions of the two partial dislocations, and 
+ to their middle, i.e. the position of the total dislocation. 

Using the same EAM potential, Khater and Bacor*2ii 
showed that a screw dislocation can also dissociate in 
the basal plane. The basal configuration is obtained by 
introducing in the same basal plane two partial disloca- 
tions of Burgers vector 1/3[1100] and 1/3[0110], in agree- 
ment with the location of the minimum on the basal 7- 
surface (Fig.[lJ. After relaxation of the atomic positions, 
the dislocation remains dissociated in the basal plane, as 
can be seen from the corresponding differential displace- 
ment map (Fig. [TU|) . A fit of the screw component of the 
disregistry created by the dislocation in the basal plane 
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leads to a dissociation length d = 6.0 A, a higher value 
than predicted by elasticity theory (Tab. [TTJ) . This con- 
figuration is metastable. It has indeed an energy higher 
by 62meVA _1 than the configuration dissociated in the 
prism plane. 

We check if ab initio calculations also lead to such a 
metastable configuration. Starting from a screw disloca- 
tion dissociated in two partial dislocations in the basal 
plane, both Siesta and Pwscf lead after relaxation of 
the atomic positions to the stable configuration dissoci- 
ated in the prism plane. This is true both with the S 
(n — 5) and the O (n — 4) periodic arrangements. These 
ab initio calculations show then that such a dissociation 
of the screw dislocation in the basal plane is unstable. 
The metastable configuration observed with the EAM 
potential appears to be an artifact of the empirical po- 
tential. This is not specific to the Mendelev and Ackland 
potential^ 4 , as a configuration dissociated in the basal 
plane is stabilized by any central forces potential ^ 5 ! 41 ' 58 

V. PEIERLS BARRIER 

Before calculating ab initio the Peierls barrier of the 
screw dislocation, we use the EAM potential to assess 
the validity of the method and to check the convergence 
of the Peierls barrier with the size of the simulation cell. 



A. Methodology 

We determine the Peierls barrier for the screw dislo- 
cation gliding in the prism plane. This is done using a 
constrained minimization between two adjacent equilib- 
rium configurations of the dislocations. We move both 
dislocations in the same direction by one Peierls distance 
Ap = c/2 between the initial and final states, so as to 
keep constant the distance between dislocations. 

Two different algorithms are used to perform the con- 
strained minimization, the simple drag method and the 
more robust nudged elastic band (NEB) method^ In- 
termediate configurations are built by linearly interpo- 
lating the atomic coordinates between the initial and fi- 
nal states. We define the corresponding reaction coor- 
dinate C = (A ~ A 1 ) • (A F - X l )/\\X F - A 1 !) 2 , where 
A, A 1 , and A F are the 3N vectors defining atomic posi- 
tions for respectively the intermediate, initial, and final 
configurations. In the drag method, the minimization is 
performed on all atomic coordinates with the constraint 
that £ remains fixed for each of the nine intermediate 
images. Nine intermediate images are also used in NEB 
method with a spring constant k = 0.1eVA~ 2 . 

So as to obtain the variation along the path of the dis- 
location energy with its position, we need to determine 
the dislocation position XDisio(C) f° r each intermediate 
image, once it has been relaxed. This is done thanks to 
a fit of the disregistry in the prism plane, as described 
in iflVBI This allows us to check that both dislocations 
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FIG. 11. Peierls barrier of a screw dislocation calculated with 
the EAM potential for the S periodic arrangement with 1600 
atoms (n = m = 20). Two different methods for finding 
the minimum energy path have been used: the NEB method 
and a simple constrained minimization (drag). Symbols cor- 
respond to the results of atomistic simulations and lines to 
their interpolation with Fourier series. 



in the simulation cell are moving in a coordinated way: 
the distance between them remains fixed. As a conse- 
quence, there is no variation of the elastic energy along 
the path and the energy variation AEp(Q obtained by 
the constrained minimization corresponds to a variation 
of the core energy, i. e. the Peierls energy. We deduce the 
Peierls energy Ai?p(a;Disio) by eliminating the reaction 
coordinate £ between A_Ep(£) and XDisio(C)- 

Fig. ITTIillustrates the whole procedure for a given dislo- 
cation periodic arrangement using both constrained min- 
imization techniques. The variation AEp(£) slightly dif- 
fers between both techniques: for a given reaction coor- 
dinate C, the drag method leads to a state of lower energy 
than NEB method. This corresponds to a small differ- 
ence in the dislocation position: this position deviates 
more with drag than with NEB method from a linear 
variation. Nevertheless, one obtains at the end the same 
Peierls barrier A£?p(a;Disio), whatever the method used. 
These differences observed for the functions AEp(() and 
a^Disio(C) between drag and NEB methods increase with 
the size of the simulation cell, i.e. with the number of 
degrees of liberty. For large simulation cells (contain- 
ing more than 3600 atoms in the S periodic arrangement 
for instance), the drag method sometimes fails to find a 
continuous path between the initial and final states: one 
has to use the NEB method in theses cases. Only much 
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and O periodic arrangements give the same Peierls bar- 
rier, and hence the same Peierls stress, for a large enough 
simulation cell (> 1000 atoms). The Peierls barrier in- 
creases when the size of the simulation cell decreases with 
the S periodic arrangement, whereas it decreases with 
the O periodic arrangement. As a consequence, the S 
and O periodic arrangement should respectively lead to 
an upper and lower limits of the Peierls stress for small 
simulation cells. We will therefore use the S periodic 
arrangement to calculate ab initio this Peierls barrier. 
This will allow us to confirm that the Peierls stress is as 
low as indicated by experiments and this EAM poten- 
tial. Moreover, it is worth pointing that the expected 
Peierls barrier should be small: AEp — 0.4meVA _1 at 
the saddle point according to the EAM potential. This 
corresponds to a difference of energy 2bAEp = 2.6 meV 
for a simulation cell of minimal height containing a dis- 
location dipole. Such a small energy difference may be 
problematic because of the precision of ab initio calcula- 
tions. Looking for an upper limit of this value is easier as 
it minimizes the problems associated with this precision. 

Finally, it is worth pointing that the dissociation length 
varies during the dislocation migration, and this variation 
is more pronounced (~ 10%) for the smallest simulation 
cells. But, like for the core energy, we did not manage to 
relate the size dependence of the Peierls barrier to this 
variation of the dissociation length. 



FIG. 12. Variation of the Peierls barrier with the size of the 
simulation cell calculated with the EAM potential for both 
periodic arrangements. 

smaller simulation cells can be studied ab initio. For 
these sizes, the drag and the NEB methods always lead 
to the same result. We will therefore only use the drag 
method in the ab initio calculations, as it costs much less 
CPU time. 

Finally, we interpolate with Fourier series the periodic 
functions AEp(() and £Disio(C) — ^pC- This leads to a 
smooth function Ai5p(a;Disio) that can be derived. The 
Peierls stress o~p is deduced from the maximal slope of 
this function, 

», fa (»**Y (9) 

b \OTDislo/ 

We obtain a Peierls stress crp = 24 ± 1 MPa for the 
EAM potential. Khater and Bacon^. determined, for the 
same empirical potential, a Peierls stress crp = 22 MPa 
using molecular statics simulations under applied stress. 
As the agreement between both methods is good, we see 
that the Peierls stress can be defined either from the slope 
of the Peierls barrier or from the minimal applied stress 
under which the dislocation glides indefinitely. 

We now examine, still with the EAM potential, how 
this Peierls barrier varies when the size of the simulation 
cell decreases up to reaching a number of atoms that can 
be handled in ab initio calculations (Fig. IT2"j) . Both the S 



B. Ab initio barriers 

The Peierls barriers obtained by ab initio calculations 
are shown on Fig. HSfa) for Pwscf and Fig. \W(b) for 
Siesta. In both cases, the height of the barrier decreases 
when the number of atoms increases, in agreement with 
what has been observed with the EAM potential for the 
same S periodic arrangement. For a given number of 
atoms, the ab initio barriers are a little bit lower than the 
ones obtained with the EAM potential. Results obtained 
with Siesta are noisy: the energy barrier that we want 
to calculate is so small that it needs a really strict conver- 
gence criterion on atomic forces for the relaxation. We 
did not manage to reach such a precision with Siesta. 
On the other hand, the barriers obtained with Pwscf 
are smooth. We can therefore interpolate the ab initio 
results with Fourier series, and estimate then the Peierls 
stress (Eq. . This leads to dp = 36 MPa for the simula- 
tion cell containing 100 atoms and crp = 21 MPa for 144 
atoms. Considering that these values, obtained in small 
simulation cells, are upper limits, ab initio calculations 
predict a Peierls stress smaller than 21 MPa for a screw 
dislocation in zirconium gliding in a prism plane. 

The comparison of this ab initio estimate of the Peierls 
stress with experimental data^ i 19 ' 57 is quite challeng- 
ing. As pointed out in the introduction, the yield stress 
of zirconium strongly depends on its oxygen content. 
No experimental data exists for a purity better than 
0.07% O (in atomic fraction). Moreover, all measure- 
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ments have been performed at temperatures higher than 
77 K. The determination of pure zirconium flow stress at 
K therefore needs some extrapolation of experimental 
data. Without a clear understanding of the mechanisms 
responsible of zirconium hardening by oxygen impurities, 
such an extrapolation is quite hazardous. Nevertheless, a 
graphical comparison fFig. H4|) of experimental data with 
our ab initio Peierls stress shows a reasonable agreement, 
keeping in mind that the value 21 MPa has to be consid- 
ered as an upper limit. 



VI. CONCLUSIONS 
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FIG. 13. Peierls barrier for a screw dislocation gliding in 
its prism plane calculated ab initio with (a) Pwscf and (b) 
Siesta for the S periodic arrangement. Symbols correspond 
to ab initio results and lines to their interpolation by Fourier 
series. 
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FIG. 14. Zirconium flow stress determined experimentally for 
various O contents.—^ ' 19 ' 57 Data have been extrapolated to 
K and are compared to the ab initio Peierls stress of screw 
dislocations in pure Zr. 



Using two ab initio approaches (Pwscf and Siesta), 
both in the DFT-GGA approximation, we have shown 
that a 1/3(1210) screw dislocation in hep zirconium dis- 
sociates in two partial dislocations with a pure screw 
character. This is in agreement with the minimum in 
1/6(1210) found for the generalized stacking fault energy 
in the prism plane. We could extract the dissociation 
length from our atomistic simulations. Although this dis- 
sociation length is small (d ^6 A for Pwscf and d ~ 4 A 
for Siesta), it is in reasonable agreement with the one 
predicted by elasticity theory. 

The EAM potential of Mendelev and Ackland^i leads 
to the same structure of the dislocation dissociated in the 
prism plane. The metastable configuration dissociated in 
the basal plane predicted by this potential is not stable 
in ab initio calculations, both with Pwscf and Siesta. 
This configuration is therefore an artifact of this poten- 
tial, probably induced by the deep minimum in 1/3(1100) 
found with this empirical potential for the generalized 
stacking fault in the basal plane. This minimum is much 
more shallow in ab initio calculations. 

We could also obtain an ab initio estimate of the 
Peierls stress of the screw dislocation gliding in the prism 
plane. Calculations with Pwscf lead to an upper limit 
of 21 MPa for this Peierls stress. This small value shows 
that screw dislocations can glide quite easily in pure zir- 
conium, thus confirming what had been obtained with 
Mendelev and Ackland EAM potential. Such a small 
Peierls stress is in agreement with experimental data, 
once the hardening of oxygen impurities has been con- 
sidered. 
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